// 试试深浅小队作品
// Created by Tong on 2020/10/25.
//

#include <iostream>
#include <windows.h> 

#include "GRDECLGrid.h"
#include "MeshInterpolation.h"



int main(int argc, char *argv[])
{   
    //获取程序当前运行的目录名称 GridA 或 GridB
    char cwd[MAX_PATH];
    std::string gridFileName;
    if (GetCurrentDirectory(MAX_PATH, cwd))
    {   
        gridFileName = cwd;
        gridFileName = gridFileName.substr(gridFileName.find_last_of("\\") + 1);
        std::cout << "Current Grid directory is: " << gridFileName << std::endl;
    }
    else
    {
        std::cerr << "Error getting Grid file name." << std::endl;
        return -1;
    }

    // 输入文件1文件名
    std::string baseGRDECLFile("baseGrid.out");
    std::cout<<"input baseGrid file:"<<baseGRDECLFile<<std::endl;
    // 输入文件1构造GRDECL原始关键字数据（解析GRDECL文件）
    GRDECLRawData rawBaseData(baseGRDECLFile);
    // 由输入文件的GRDECL关键字数据构造网格
    GRDECLGrid grdeclBaseGrid(rawBaseData);
    // 输出输入文件的vtk可视化文件
    grdeclBaseGrid.vtkCellWriter("baseGrid.vtk");
    IntList baseLayer;

    // 输出断层面几何数据
    //grdeclBaseGrid.vtkFaultFaceWriter("baseGRDECLFaultFace.vtk");
    // // 构建延拓与剖分约束几何
    // grdeclBaseGrid.genBoundaryFace();
    // // 输出延拓与剖分约束几何
    // grdeclBaseGrid.vtkBoundaryWriter("baseGRDECLBoundary2.vtk");
    //grdeclBaseGrid.freeData2();

    // 输入文件2文件名
    std::string fineGRDECLFile("fineGrid.out");
    std::cout<<"input fineGrid file:"<<fineGRDECLFile<<std::endl;
    GRDECLRawData rawFineData(fineGRDECLFile,"Props.txt");
    // 由输入文件的GRDECL关键字数据构造网格
    GRDECLGrid::extendZMax_ = 1000;
    GRDECLGrid::extendZMin_ = -1200;
    if (gridFileName == "GridA")
    {
        GRDECLGrid::extendZMax_ = 1000;
        GRDECLGrid::extendZMin_ = -1200;
        baseLayer={10,14,15,21};
    }
    else if (gridFileName == "GridB")
    {
        GRDECLGrid::extendZMax_ = -400;
        GRDECLGrid::extendZMin_ = -1050;
        baseLayer={10,17,18,20,32};
    }
    else
    {
        std::cerr << "Unknown Grid file name!" << std::endl;
        return -1;
    }


    GRDECLGrid grdeclFineGrid(rawFineData);
    // GRDECLGrid 守恒属性数据
    const double totPoroVolume=GRDECLGrid::calcConservation(grdeclFineGrid.Poro_,grdeclFineGrid.cellVolumes_,grdeclFineGrid.Poro_.size());
    const double totMass = GRDECLGrid::calcConservation(grdeclFineGrid.RockDens_,grdeclFineGrid.cellVolumes_,grdeclFineGrid.RockDens_.size());
    std::cout<<"Total Poro Volume: "<<totPoroVolume<<std::endl;
    std::cout<<"Total Mass: "<<totMass<<std::endl;
    // 输出输入文件的vtk可视化文件
    grdeclFineGrid.vtkCellWriter("fineGrid.vtk");
    // 输出断层面几何数据
    //grdeclFineGrid.vtkFaultFaceWriter("fineGRDECLFaultFace.vtk");
    // 边界输出
    //grdeclFineGrid.vtkBoundaryWriter("fineGRDECLBoundary.vtk");
    // 六面体网格分解
    std::cout<<"--------Start decompose Hex Cell--------"<<std::endl;
    clock_t startT = std::clock();
    grdeclFineGrid.decomposeAllHexCell();
    //grdeclFineGrid.vtkHexCellWriter("fineGridHex.vtk");
    //grdeclFineGrid.vtkKMaxWriter("fineGRDECLKMax.vtk");
    // 延拓六面体网格
    grdeclFineGrid.genHexExtendCell();
    std::cout<< "Decompose and extend hex cell using: "<<double(std::clock()-startT)/CLOCKS_PER_SEC<<"s"<<std::endl;
    // grdeclFineGrid.vtkExtendedCell(
    //     "fineGridHexExtend.vtk",
    //     grdeclFineGrid.allHexCoord_,
    //     grdeclFineGrid.allHexCellNodes_,
    //     grdeclFineGrid.hexPointData_,
    //     grdeclFineGrid.hexCellData_);
    // 网格插值与属性设置
    // 构造插值器 
    clock_t interStartT = std::clock();
    const size_t hexCellNum = grdeclFineGrid.hexCellCenter_.size();
    const size_t hexReserviorCellNum = int(grdeclFineGrid.hexReserviorCellNum_);
    MeshInterpolation interpolatorHex(
        grdeclFineGrid.cellCenters_,
        grdeclFineGrid.hexCellCenter_,
        grdeclFineGrid.cellCenters_.size(),
        hexReserviorCellNum);
    // 油藏区域插值网格1到网格2
    interpolatorHex.interpolateMesh1toMesh2(grdeclFineGrid.Poro_, grdeclFineGrid.hexPoro_,hexReserviorCellNum);
    interpolatorHex.interpolateMesh1toMesh2(grdeclFineGrid.RockDens_, grdeclFineGrid.hexRockDens_,hexReserviorCellNum);
    interpolatorHex.interpolateMesh1toMesh2(grdeclFineGrid.YoungsMod_, grdeclFineGrid.hexYoungsMod_,hexReserviorCellNum);
    interpolatorHex.interpolateMesh1toMesh2(grdeclFineGrid.Poisson_, grdeclFineGrid.hexPoisson_,hexReserviorCellNum);
    std::cout<< "Interpolation hex data using: "<<double(std::clock()-interStartT)/CLOCKS_PER_SEC<<"s"<<std::endl;

    // 延拓区属性赋值
    // 赋值延拓区孔隙度
    GRDECLGrid::assignCellData(
        grdeclFineGrid.hexCellData_,
        grdeclFineGrid.hexPoro_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区岩石密度
    GRDECLGrid::assignCellData(
        grdeclFineGrid.hexCellData_,
        grdeclFineGrid.hexPoro_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区杨氏模量
    GRDECLGrid::assignCellData(
        grdeclFineGrid.hexCellData_,
        grdeclFineGrid.hexYoungsMod_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区泊松比
    GRDECLGrid::assignCellData(
        grdeclFineGrid.hexCellData_,
        grdeclFineGrid.hexPoisson_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 提取地层面
    GRDECLGrid::extractLayer(grdeclFineGrid.hexPointData_,baseLayer);

    const double hexTotPoroVolume = GRDECLGrid::calcConservation(grdeclFineGrid.hexPoro_, grdeclFineGrid.hexCellVolume_, hexReserviorCellNum);
    const double hexTotalMass = GRDECLGrid::calcConservation(grdeclFineGrid.hexRockDens_, grdeclFineGrid.hexCellVolume_, hexReserviorCellNum);

    std::cout << "Hex Total PoroVolume Scale Factor: " << totPoroVolume / hexTotPoroVolume << std::endl;
    std::cout << "Hex Total Mass Scale Factor: " << totMass / hexTotalMass << std::endl;

    GRDECLGrid::scaleData(grdeclFineGrid.hexPoro_, hexReserviorCellNum, totPoroVolume / hexTotPoroVolume);
    GRDECLGrid::scaleData(grdeclFineGrid.hexRockDens_, hexReserviorCellNum, totMass / hexTotalMass);

    //输出六面体分解的所有数据
    GRDECLGrid::vtkExtendedCellData(
        "fineGridHexExtendData.vtk",
        grdeclFineGrid.allHexCoord_,
        grdeclFineGrid.allHexCellNodes_,
        grdeclFineGrid.hexPointData_,
        grdeclFineGrid.hexCellData_,
        grdeclFineGrid.hexYoungsMod_,
        grdeclFineGrid.hexPoisson_,
        grdeclFineGrid.hexPoro_,
        grdeclFineGrid.hexRockDens_);
    
    // 释放Hex网格数据
    //system("pause for memory check");
    grdeclFineGrid.freeHexData();
    interpolatorHex.freeAll();
    std::cout << "--------End decompose Hex Cell Using " << double(std::clock() - startT) / CLOCKS_PER_SEC << "s --------" << std::endl;

    // // 全四面体网格分解
    std::cout << "--------Start decompose Tet Cell--------" << std::endl;
    startT = std::clock();
    grdeclFineGrid.decomposeAllTetCell();
    grdeclFineGrid.genTetExtendCell();
    std::cout << "Decompose and extend tet cell using: " << double(std::clock() - startT) / CLOCKS_PER_SEC << "s" << std::endl;
    // grdeclFineGrid.vtkExtendedCell(
    //     "fineGridTetExtend.vtk",
    //     grdeclFineGrid.allTetCoord_,
    //     grdeclFineGrid.allTetCellNodes_,
    //     grdeclFineGrid.tetPointData_,
    //     grdeclFineGrid.tetCellData_);
    // 网格插值与属性设置
    // 构造插值器
    interStartT = std::clock();
    const size_t tetCellNum = grdeclFineGrid.tetCellCenter_.size();
    const size_t tetReserviorCellNum = int(grdeclFineGrid.tetReserviorCellNum_);
    MeshInterpolation interpolatorTet(
        grdeclFineGrid.cellCenters_,
        grdeclFineGrid.tetCellCenter_,
        grdeclFineGrid.cellCenters_.size(),
        tetReserviorCellNum);
    // 油藏区域插值网格1到网格2
    interpolatorTet.interpolateMesh1toMesh2(grdeclFineGrid.Poro_, grdeclFineGrid.tetPoro_,tetReserviorCellNum);
    interpolatorTet.interpolateMesh1toMesh2(grdeclFineGrid.RockDens_, grdeclFineGrid.tetRockDens_,tetReserviorCellNum);
    interpolatorTet.interpolateMesh1toMesh2(grdeclFineGrid.YoungsMod_, grdeclFineGrid.tetYoungsMod_,tetReserviorCellNum);
    interpolatorTet.interpolateMesh1toMesh2(grdeclFineGrid.Poisson_, grdeclFineGrid.tetPoisson_,tetReserviorCellNum);
    std::cout << "Interpolation tet data using: " << double(std::clock() - interStartT) / CLOCKS_PER_SEC << "s" << std::endl;
    // 延拓区属性赋值
    // 赋值延拓区孔隙度
    GRDECLGrid::assignCellData(
        grdeclFineGrid.tetCellData_,
        grdeclFineGrid.tetPoro_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区岩石密度
    GRDECLGrid::assignCellData(
        grdeclFineGrid.tetCellData_,
        grdeclFineGrid.tetPoro_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区杨氏模量
    GRDECLGrid::assignCellData(
        grdeclFineGrid.tetCellData_,
        grdeclFineGrid.tetYoungsMod_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区泊松比
    GRDECLGrid::assignCellData(
        grdeclFineGrid.tetCellData_,
        grdeclFineGrid.tetPoisson_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 提取地层面
    GRDECLGrid::extractLayer(grdeclFineGrid.tetPointData_,baseLayer);
    
    
    const double tetTotPoroVolume = GRDECLGrid::calcConservation(grdeclFineGrid.tetPoro_, grdeclFineGrid.tetCellVolume_, tetReserviorCellNum);
    const double tetTotalMass = GRDECLGrid::calcConservation(grdeclFineGrid.tetRockDens_, grdeclFineGrid.tetCellVolume_, tetReserviorCellNum);

    std::cout << "Tet Total PoroVolume Scale Factor:" << totPoroVolume / tetTotPoroVolume << std::endl;
    std::cout << "Tet Total Mass Scale Factor:" << totMass / tetTotalMass << std::endl;
    GRDECLGrid::scaleData(grdeclFineGrid.tetPoro_, tetReserviorCellNum, totPoroVolume / tetTotPoroVolume);
    GRDECLGrid::scaleData(grdeclFineGrid.tetRockDens_, tetReserviorCellNum, totMass / tetTotalMass);

    //输出四面体分解的所有数据
    GRDECLGrid::vtkExtendedCellData(
        "fineGridTetExtendData.vtk",
        grdeclFineGrid.allTetCoord_,
        grdeclFineGrid.allTetCellNodes_,
        grdeclFineGrid.tetPointData_,
        grdeclFineGrid.tetCellData_,
        grdeclFineGrid.tetYoungsMod_,
        grdeclFineGrid.tetPoisson_,
        grdeclFineGrid.tetPoro_,
        grdeclFineGrid.tetRockDens_);

    // 释放Tet网格数据
    //system("pause for memory check");
    grdeclFineGrid.freeTetData();
    interpolatorTet.freeAll();
    std::cout << "--------End decompose Tet Cell Using " << double(std::clock() - startT) / CLOCKS_PER_SEC << "s --------" << std::endl;


    
    // grdeclFineGrid.vtkTetCellWriter("fineGridTet.vtk");
    // grdeclFineGrid.decomposeCell();
    // grdeclFineGrid.vtkDecomposeCellWriter("fineGridDecompose.vtk");
    
    
    // // 全六面体分解网格
    // grdeclGrid.decomposeAllHexCell();
    // // 输出全六面体分解后网格几何
    // grdeclGrid.vtkHexCellWriter("GRDECLHexGrid.vtk");
    // system("pause for memory check");
    // // 全四面体分解网格
    // grdeclGrid.decomposeAllTetCell();
    // // 输出全四面体分解后网格几何
    // grdeclGrid.vtkTetCellWriter("GRDECLTetGrid.vtk");
    // system("pause for memory check");
    // 四面体-金字塔-六面体混合网格分解
    std::cout << "--------Start decompose Hybrid Cell--------" << std::endl;
    startT = std::clock();
    grdeclFineGrid.decomposeCell();
    // 延拓混合网格
    grdeclFineGrid.genHybridExtendCell();
    std::cout << "Decompose and extend hybrid cell using: " << double(std::clock() - startT) / CLOCKS_PER_SEC << "s" << std::endl;

    // grdeclFineGrid.vtkExtendedCell(
    //     "fineGridHybridExtend.vtk",
    //     grdeclFineGrid.decomposedCoord_,
    //     grdeclFineGrid.decomposedCellNodes_,
    //     grdeclFineGrid.decomposedPointData_,
    //     grdeclFineGrid.decomposedCellData_);
    // 网格插值与属性设置
    // 构造插值器 
    interStartT = std::clock();
    const size_t hybridCellNum = grdeclFineGrid.decomposedCellCenter_.size();
    const size_t hybridReserviorCellNum = int(grdeclFineGrid.decompReservoirCellNum_);
    MeshInterpolation interpolatorHybrid(
        grdeclFineGrid.cellCenters_,
        grdeclFineGrid.decomposedCellCenter_,
        grdeclFineGrid.cellCenters_.size(),
        hybridReserviorCellNum);
    // // 输出混合分解后网格几何
    // 油藏区域插值网格1到网格2
    interpolatorHybrid.interpolateMesh1toMesh2(grdeclFineGrid.Poro_, grdeclFineGrid.decomposedPoro_,hybridReserviorCellNum);
    interpolatorHybrid.interpolateMesh1toMesh2(grdeclFineGrid.RockDens_, grdeclFineGrid.decomposedRockDens_,hybridReserviorCellNum);
    interpolatorHybrid.interpolateMesh1toMesh2(grdeclFineGrid.YoungsMod_, grdeclFineGrid.decomposedYoungsMod_,hybridReserviorCellNum);
    interpolatorHybrid.interpolateMesh1toMesh2(grdeclFineGrid.Poisson_, grdeclFineGrid.decomposedPoisson_,hybridReserviorCellNum);
    std::cout << "Interpolation hybrid data using: " << double(std::clock() - interStartT) / CLOCKS_PER_SEC << "s" << std::endl;
    // 延拓区属性赋值
    // 赋值延拓区孔隙度
    GRDECLGrid::assignCellData(
        grdeclFineGrid.decomposedCellData_,
        grdeclFineGrid.decomposedPoro_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区岩石密度
    GRDECLGrid::assignCellData(
        grdeclFineGrid.decomposedCellData_,
        grdeclFineGrid.decomposedPoro_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区杨氏模量
    GRDECLGrid::assignCellData(
        grdeclFineGrid.decomposedCellData_,
        grdeclFineGrid.decomposedYoungsMod_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 赋值延拓区泊松比
    GRDECLGrid::assignCellData(
        grdeclFineGrid.decomposedCellData_,
        grdeclFineGrid.decomposedPoisson_,
        0.1,  // over
        0.2,  // under
        0.3,  // side
        0.4,  // overside
        0.5); // underside
    // 提取地层面
    GRDECLGrid::extractLayer(grdeclFineGrid.decomposedPointData_,baseLayer);

    const double hybridTotPoroVolume = GRDECLGrid::calcConservation(grdeclFineGrid.decomposedPoro_, grdeclFineGrid.decomposedCellVolume_, hybridReserviorCellNum);
    const double hybridTotalMass = GRDECLGrid::calcConservation(grdeclFineGrid.decomposedRockDens_, grdeclFineGrid.decomposedCellVolume_, hybridReserviorCellNum);

    std::cout << "Hybrid Total PoroVolume Scale Factor:" << totPoroVolume / hybridTotPoroVolume << std::endl;
    std::cout << "Hybrid Total Mass Scale Factor:" << totMass / hybridTotalMass << std::endl;
    GRDECLGrid::scaleData(grdeclFineGrid.decomposedPoro_, hybridReserviorCellNum, totPoroVolume / hybridTotPoroVolume);
    GRDECLGrid::scaleData(grdeclFineGrid.decomposedRockDens_, hybridReserviorCellNum, totMass / hybridTotalMass);

    //输出混合分解的所有数据
    GRDECLGrid::vtkExtendedCellData(
        "fineGridHybridExtendData.vtk",
        grdeclFineGrid.decomposedCoord_,
        grdeclFineGrid.decomposedCellNodes_,
        grdeclFineGrid.decomposedPointData_,
        grdeclFineGrid.decomposedCellData_,
        grdeclFineGrid.decomposedYoungsMod_,
        grdeclFineGrid.decomposedPoisson_,
        grdeclFineGrid.decomposedPoro_,
        grdeclFineGrid.decomposedRockDens_);
    //system("pause for memory check");
    std::cout << "--------End decompose Hybrid Cell Using " << double(std::clock() - startT) / CLOCKS_PER_SEC << "s --------" << std::endl;
    return 0;
}
